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HIGH  ORDER  FINITE  DIFFERENCE  AND  FINITE  VOLUME  WENO  SCHEMES  AND 
DISCONTINUOUS  GALERKIN  METHODS  FOR  CFD 

CHI- WANG  SHU* 

Abstract.  In  recent  years  high  order  numerical  methods  have  been  widely  used  in  computational  fluid 
dynamics  (CFD),  to  effectively  resolve  complex  flow  features  using  meshes  which  are  reasonable  for  today’s 
computers.  In  this  paper  we  review  and  compare  three  types  of  high  order  methods  being  used  in  CFD, 
namely  the  weighted  essentially  non-oscillatory  (WENO)  finite  difference  methods,  the  WENO  finite  volume 
methods,  and  the  discontinuous  Galerkin  (DG)  finite  element  methods.  We  summarize  the  main  features 
of  these  methods,  from  a  practical  user’s  point  of  view,  indicate  their  applicability  and  relative  strength, 
and  show  a  few  selected  numerical  examples  to  demonstrate  their  performance  on  illustrative  model  CFD 
problems. 

Key  words,  weighted  essentially  non-oscillatory,  discontinuous  Galerkin,  finite  difference  method,  finite 
volume  method,  computational  fluid  dynamics 

Subject  classification.  Applied  and  Numerical  Mathematics 

1 .  Introduction.  In  recent  years  high  order  numerical  methods  have  been  widely  used  in  computational 
fluid  dynamics  (CFD),  to  effectively  resolve  complex  flow  features.  In  this  paper  we  refer  to  high  order 
methods  by  those  with  order  of  accuracy  at  least  three.  Traditionally,  first  and  second  order  numerical 
methods  are  often  preferred  in  practical  calculations,  because  of  their  simplicity  and  robustness  (i.e.  one  can 
always  get  some  output,  although  it  may  not  be  very  accurate).  On  the  other  hand,  high  order  methods 
often  give  the  impression  of  being  complicated  to  understand  and  to  code,  and  costly  to  run  (on  the  same 
mesh  compared  with  lower  order  methods),  and  less  robust  (the  code  may  blow  up  in  tough  situations  when 
the  lower  order  methods  still  give  stable  output).  In  this  paper,  we  hope  to  at  least  partially  dispel  this 
impression  about  high  order  methods,  using  three  typical  types  of  high  order  methods  as  examples. 

Before  we  move  on  to  the  details  of  high  order  methods,  let  us  point  out  that,  at  least  in  certain 
situations,  the  solution  structures  are  so  complicated  and  the  time  of  evolution  of  these  structures  so  long 
that  it  is  impractical  to  use  low  order  methods  to  obtain  an  acceptable  resolution.  Often  such  problems 
involve  both  shocks  and  complicated  smooth  region  structures,  calling  for  special  non-oscillatory  type  high 
order  schemes  which  are  emphasized  in  this  paper.  A  very  simple  example  to  illustrate  this  is  the  evolution 
of  a  two  dimensional  periodic  vortex  for  the  compressible  Euler  equations,  which  was  first  used  in  [41],  One 
could  add  a  shock  to  this  problem  so  that  it  becomes  a  problem  of  shock  interaction  with  a  vortex,  which 
is  very  typical  in  aeroacoustics.  However  we  will  consider  the  solution  without  a  shock  here  as  it  admits 
an  analytically  given  exact  solution,  making  it  easier  to  compare  numerical  resolutions  of  different  schemes. 
The  Euler  equations  are  given  by  a  conservation  law 

ut  +  f(u)x  +  g(u)y  =  0,  (1-1.1) 


‘Division  of  Applied  Mathematics,  Brown  University,  Providence,  RI  02912.  E-mail:  shu@cfm.brown.edu.  Research  sup¬ 
ported  by  ARO  grants  DAAG55-97-1-0318  and  DAAD19-00-1-0405,  NSF  grants  DMS-9804985  and  ECS-9906606,  NASA  Lang¬ 
ley  grant  NAG-1-2070  and  Contract  NAS1-97046  while  the  author  was  in  residence  at  ICASE,  NASA  Langley  Research  Center, 
Hampton,  VA  23681-2199,  and  AFOSR  grant  F49620-99- 1-0077. 


1 


where 


u  =  (p,pu,i,pu2,E), 


f(u)  =  ( pu i ,  pu\  +  p,  puiu-2  ,Ui(E  +p)),  g(u )  =  (pu2 ,  puiu2 ,  pu\  +  p,  u2(E  +  p)), 

Here  p  is  the  density,  {u\,u2)  is  the  velocity,  E  is  the  total  energy,  p  is  the  pressure,  related  to  the  total 
energy  E  by 

E  =  +  ~p{u\  +  Wo) 

7  —  1  2 

with  7  =  1.4  for  air.  The  periodic  vortex  problem  is  set  up  in  a  computational  domain  [0,10] x [0,10].  The 
boundary  condition  is  periodic  in  both  directions.  The  initial  condition  is  given  by 

(u1{x,y,0),u,2{x,y,0))  =  (1,1)  +  ^-e°-5{1~r2)(-y,x), 


0)  =  1  - 


(7  ~  l)f2  ,1-r2 
877T2 


S(x,y,  0)  =  1, 


where  the  temperature  T  and  the  entropy  S  are  related  to  the  density  p  and  the  pressure  p  by 


S  = 


P_ 

P 7 


and  (x,y)  =  { x  —  5,y  —  5),  r2  =  x2+y2,  and  the  vortex  strength  e  =  5. 

It  can  be  readily  verified  that  the  Euler  equations  with  the  above  initial  conditions  admit  an  exact 
solution  which  is  converted  with  the  speed  (1,1)  in  the  diagonal  direction.  Because  of  the  periodic  boundary 
condition,  we  can  simulate  this  flow  for  a  very  long  time.  We  first  show  the  simulation  results  at  t  =  10, 
namely  after  one  time  period.  When  we  perform  the  simulation  with  a  second  order  finite  difference  MUSCL 
type  TVD  scheme  and  a  fifth  order  finite  difference  WENO  scheme  with  the  same  uniform  mesh  of  802  points, 
for  this  relatively  short  time,  although  the  second  order  scheme  gives  inferior  results  comparing  with  that 
of  the  fifth  order  scheme,  Fig.  1.1,  it  may  be  argued  that  the  second  order  scheme  still  gives  an  acceptable 
resolution.  If  we  increase  the  number  of  mesh  points  for  the  second  order  scheme  to  2002  points,  see  Fig.  1.3, 
left,  then  the  resolution  is  roughly  comparable  to  that  of  the  fifth  order  WENO  scheme  using  802  points  in 
Fig.  1.1,  right.  A  two  dimensional  time  dependent  simulation  with  a  2002  mesh  has  2.53  =  15.6  times  more 
space  time  mesh  points  than  a  802  mesh.  Considering  that  the  CPU  time  of  a  fifth  order  finite  difference 
WENO  scheme  is  roughly  3  to  8  times  more  than  that  of  a  second  order  TVD  scheme  on  the  same  mesh 
(depending  on  the  specific  forms  of  the  schemes  and  time  discretization),  we  could  conclude  that  the  second 
order  TVD  scheme  has  a  larger  but  still  comparable  CPU  cost  than  the  fifth  order  WENO  scheme  to  reach 
the  same  resolution,  for  this  problem  with  relatively  short  time.  When  we  look  at  the  result  at  t  =  100, 
namely  after  ten  time  periods,  the  situation  changes  dramatically.  On  the  same  802  mesh,  one  can  see  in 
Fig.  1.2  that  the  second  order  finite  difference  TVD  scheme  has  a  much  worse  resolution  than  the  fifth  order 
finite  difference  WENO  scheme.  Clearly  the  result  of  the  second  order  scheme  with  this  mesh  for  this  long 
time  is  completely  unacceptable.  This  time,  even  if  one  increases  the  number  of  mesh  points  to  3202  for  the 
second  order  scheme  (which  makes  the  CPU  time  for  such  a  run  magnitudes  more  than  that  of  a  fifth  order 
WENO  scheme  on  a  802  mesh),  it  still  does  not  provide  a  satisfactory  resolution,  Fig.  1.3,  right.  A  more 
refined  mesh  would  not  be  practical  for  three  dimensional  simulations.  Clearly,  in  this  situation  the  second 
order  scheme  is  inadequate  to  provide  a  satisfactory  resolution  within  the  limit  of  today’s  computer. 
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Fig.  1.1.  Vortex  evolution.  Cut  at  x  =  5.  Density  p.  802  uniform  mesh,  t  =  10  (after  one  time  period).  Solid:  exact 
solution ;  circles:  computed  solution.  Left:  second  order  TVD  scheme:  right:  fifth  order  WENO  scheme. 


Fig.  1.2.  Vortex  evolution.  Cut  at  x  =  5.  Density  p.  802  uniform  mesh,  t  =  100  (after  10  time  periods).  Solid:  exact 
solution;  circles:  computed  solution.  Left:  second  order  TVD  scheme ;  right:  fifth  order  WENO  scheme. 

There  are  many  high  order  methods  being  used  in  CFD.  In  this  paper  we  only  discuss  three  types  of 
them: 

1.  The  weighted  essentially  non-oscillatory  (WENO)  finite  difference  methods; 

2.  The  WENO  finite  volume  methods; 

3.  The  discontinuous  Galerkin  (DG)  finite  element  methods. 

These  are  methods  suitable  for  solving  hyperbolic  conservation  laws,  such  as  the  compressible  Euler  equations 
(1.1.1),  or  convection  dominated  convection  diffusion  problems,  such  as  the  compressible  Navier-Stokes 
equations  with  high  Reynolds  numbers.  For  such  problems  shocks  and  other  discontinuities  or  high  gradient 
regions  exist  in  the  solutions,  making  it  difficult  to  design  stable  and  high  order  numerical  methods. 
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Fig.  1.3.  Vortex  evolution.  Cut  at  x  =  5.  Density  p.  Second  order  TVD  scheme.  Solid:  exact  solution ;  circles:  computed 
solution.  Left:  2002  uniform  mesh,  t  =  10  (after  one  time  period);  right:  3202  uniform  mesh,  t  =  100  (after  10  time  periods). 

Let  us  first  give  some  historical  remarks  about  these  methods. 

WENO  schemes  are  designed  based  on  the  successful  essentially  non-oscillatory  (ENO)  schemes  in  [19, 
43,  44].  The  first  WENO  scheme  is  constructed  in  [31]  for  a  third  order  finite  volume  version  in  one  space 
dimension.  In  [23],  third  and  fifth  order  finite  difference  WENO  schemes  in  multi  space  dimensions  are 
constructed,  with  a  general  framework  for  the  design  of  the  smoothness  indicators  and  nonlinear  weights. 
Later,  second,  third  and  fourth  order  finite  volume  WENO  schemes  for  2D  general  triangulation  have  been 
developed  in  [14]  and  [20].  Very  high  order  finite  difference  WENO  schemes  (for  orders  between  7  and  11) 
have  been  developed  in  [2].  Central  WENO  schemes  have  been  developed  in  [25],  [26]  and  [27]. 

Both  ENO  and  WENO  use  the  idea  of  adaptive  stencils  in  the  reconstruction  procedure  based  on  the 
local  smoothness  of  the  numerical  solution  to  automatically  achieve  high  order  accuracy  and  non-oscillatory 
property  near  discontinuities.  ENO  uses  just  one  (optimal  in  some  sense)  out  of  many  candidate  stencils 
when  doing  the  reconstruction;  while  WENO  uses  a  convex  combination  of  all  the  candidate  stencils,  each 
being  assigned  a  nonlinear  weight  which  depends  on  the  local  smoothness  of  the  numerical  solution  based 
on  that  stencil.  WENO  improves  upon  ENO  in  robustness,  better  smoothness  of  fluxes,  better  steady  state 
convergence,  better  provable  convergence  properties,  and  more  efficiency.  For  more  details  of  ENO  and 
WENO  schemes,  we  refer  to  the  lecture  notes  [41,  42], 

WENO  schemes  have  been  widely  used  in  applications.  Some  of  the  examples  include  dynamical  response 
of  a  stellar  atmosphere  to  pressure  perturbations  [13];  shock  vortex  interactions  and  other  gas  dynamics  prob¬ 
lems  [17],  [18];  incompressible  flow  problems  [47];  Hamilton-Jacobi  equations  [21];  magneto-hydrodynamics 
[24] ;  underwater  blast-wave  focusing  [28] ;  the  composite  schemes  and  shallow  water  equations  [29] ,  [30] ,  real 
gas  computations  [32],  wave  propagation  using  Fey’s  method  of  transport  [33];  etc. 

Discontinuous  Galerkin  (DG)  methods  are  a  class  of  finite  element  methods  using  completely  discontin¬ 
uous  basis  functions,  which  are  usually  chosen  as  piecewise  polynomials.  Since  the  basis  functions  can  be 
completely  discontinuous,  these  methods  have  the  flexibility  which  is  not  shared  by  typical  finite  element 
methods,  such  as  the  allowance  of  arbitrary  triangulation  with  hanging  nodes,  complete  freedom  in  changing 
the  polynomial  degrees  in  each  element  independent  of  that  in  the  neighbors  (p  adaptivity),  and  extremely 
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local  data  structure  (elements  only  communicate  with  immediate  neighbors  regardless  of  the  order  of  ac¬ 
curacy  of  the  scheme)  and  the  resulting  embarrassingly  high  parallel  efficiency  (usually  more  than  99%  for 
a  fixed  mesh,  and  more  than  80%  for  a  dynamic  load  balancing  with  adaptive  meshes  which  change  often 
during  time  evolution),  see,  e.g.  [4],  A  very  good  example  to  illustrate  the  capability  of  the  discontinuous 
Galerkin  method  in  h-p  adaptivity,  efficiency  in  parallel  dynamic  load  balancing,  and  excellent  resolution 
properties  is  the  recent  successful  simulation  of  the  Rayleigh- Taylor  flow  instabilities  in  [36]. 

The  first  discontinuous  Galerkin  method  was  introduced  in  1973  by  Reed  and  Hill  [35],  in  the  frame¬ 
work  of  neutron  transport,  i.e.  a  time  independent  linear  hyperbolic  equation.  A  major  development  of  the 
DG  method  is  carried  out  by  Cockburn  et  al.  in  a  series  of  papers  [6,  7,  8,  9],  in  which  they  have  estab¬ 
lished  a  framework  to  easily  solve  nonlinear  time  dependent  problems,  such  as  the  Euler  equations  (1.1.1), 
using  explicit,  nonlinearly  stable  high  order  Runge-Kutta  time  discretizations  [43]  (see  section  2)  and  DG 
discretization  in  space  with  exact  or  approximate  Riemann  solvers  as  interface  fluxes  and  total  variation 
bounded  (TVB)  nonlinear  limiters  to  achieve  non-oscillatory  properties  for  strong  shocks. 

The  DG  method  has  found  rapid  applications  in  such  diverse  areas  as  aeroacoustics,  electro-magnetism, 
gas  dynamics,  granular  flows,  magneto-hydrodynamics,  meteorology,  modeling  of  shallow  water,  oceanogra¬ 
phy,  oil  recovery  simulation,  semiconductor  device  simulation,  transport  of  contaminant  in  porous  media, 
turbomachinery,  turbulent  flows,  viscoelastic  flows  and  weather  forecasting,  among  many  others.  For  more 
details,  we  refer  to  the  survey  paper  [12],  and  other  papers  in  that  Springer  volume,  which  contains  the 
conference  proceedings  of  the  First  International  Symposium  on  Discontinuous  Galerkin  Methods  held  at 
Newport,  Rhode  Island  in  1999.  The  extensive  review  paper  [11]  is  also  a  good  reference  for  many  details. 

This  paper  is  written  from  a  practical  user’s  point  of  view.  We  will  not  emphasize  the  discussion  of 
theoretical  properties  of  the  schemes.  Rather,  we  will  indicate  the  practical  aspects  in  the  implementation 
of  the  algorithms,  their  applicability  in  different  situations,  and  their  relative  advantages,  and  present  a  few 
selected  numerical  examples  to  demonstrate  their  performance  on  illustrative  model  CFD  problems. 

2.  Time  discretizations.  Before  discussing  the  spatial  discretizations,  let  us  first  discuss  the  time 
discretization.  For  all  three  types  of  spatial  discretizations  discussed  in  this  paper,  we  shall  use  the  same 
time  discretization,  namely  a  class  of  high  order  nonlinearly  stable  Runge-Kutta  time  discretizations.  A 
distinctive  feature  of  this  class  of  time  discretizations  is  that  they  are  convex  combinations  of  first  order 
forward  Euler  steps,  hence  they  maintain  strong  stability  properties  in  any  semi-norm  (total  variation  norm, 
maximum  norm,  entropy  condition,  etc.)  of  the  forward  Euler  step.  Thus  one  only  needs  to  prove  nonlinear 
stability  for  the  first  order  forward  Euler  step,  which  is  relatively  easy  in  many  situations  (e.g.  TVD  schemes), 
and  one  automatically  obtains  the  same  strong  stability  property  for  the  higher  order  time  discretizations 
in  this  class.  These  methods  were  first  developed  in  [43]  and  [40],  and  later  generalized  in  [15]  and  [16].  The 
most  popular  scheme  in  this  class  is  the  following  third  order  Runge-Kutta  method  for  solving 

ut  =  L{u,t) 

where  L(u,t )  is  a  spatial  discretization  operator  (it  does  not  need  to  be,  and  often  is  not,  linear!): 

u{1)  =un  +AtL(un,tn) 

u®]  =  +  7AiL(u(1),in  +  At) 

4  4  4 

(/',j  1  =  \un  +  |w(2)  +  |  A  tL(u{2),tn  +  l  At). 

o  o  o  £ 

All  the  numerical  examples  presented  in  this  paper  are  obtained  with  this  Runge-Kutta  time  discretization. 
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3.  Finite  difference  WENO  schemes.  A  conservative  finite  difference  spatial  discretization  to  a 
conservation  law  such  as  (1.1.1)  approximates  the  derivative  f(u)x  by  a  conservative  difference 

f(u)x\x=xj  ~  (jj+ 1/2  ~  fj- 1/2  j 

where  fj+1/2  is  the  numerical  flux,  which  typically  is  a  Lipschitz  continuous  function  of  several  neighboring 
values  Ui .  g(u)y  is  approximated  in  the  same  way.  Hence  finite  difference  methods  have  the  same  format 
for  one  and  several  space  dimensions,  which  is  a  big  advantage.  For  the  simplest  case  of  a  scalar  equation 
(1.1.1)  and  if  f'(u)  >  0,  the  fifth  order  finite  difference  WENO  scheme  has  the  flux  given  by 

fj+l/2  =  «h/j+1/2  +  ^o/jji/o  +  Wj+J/o 

where  are  three  third  order  fluxes  on  three  different  stencils  given  by 

fj+1/2  =  2)  - \f^uJ- 1)  + 

/]+ 1/2  =  1)  +  \f(uj)  +  \f(uj+ 1), 

fj+1/2  =  +  ^f(uj+ 1)  -  lf(Uj+2), 

and  the  nonlinear  weights  w,;  are  given  by 


Bi  =  (f(uj- 2)  ~  2/(uj_i)  +  f(uj))2  +  ^  (f(uj- 2)  -  !/(«,  1)  +  3 f(uj))2 
If  =  (f(uj- 1)  -  2 f(Uj)  +  f{uj+ 1))2  +  i  (f(uj- 1)  -  f(uj+ 1))2 

h  |t7  (/(W'i)  -  2/(wj+i)  +  /(Wj+2))2  +  ^  (3/(uj)  -  4/(mj+i)  +  f(uJ+2))2  ■ 

Finally,  e  is  a  parameter  to  avoid  the  denominator  to  become  0  and  is  usually  taken  as  e  =  10-6  in  the 
computation. 

This  finishes  the  description  of  the  fifth  order  finite  difference  WENO  scheme  [23]  in  the  simplest  case. 
As  we  can  see,  the  algorithm  is  actually  quite  simple  and  there  are  no  tunable  parameters  in  the  scheme. 

We  summarize  the  properties  of  this  WENO  finite  difference  scheme.  For  details  of  proofs  and  numerical 
verifications,  see  [23]  and  [41,  42], 

1.  The  scheme  is  proven  to  be  uniformly  fifth  order  accurate  including  at  smooth  extrema,  and  this  is 
verified  numerically. 

2.  Near  discontinuities  the  scheme  behaves  very  similarly  to  an  ENO  scheme  [19,  43,  44],  namely  the 
solution  has  a  sharp  and  non-oscillatory  discontinuity  transition. 

3.  The  numerical  flux  has  the  same  smoothness  dependency  on  its  arguments  as  that  of  the  physical 
flux  f(u).  This  helps  in  a  convergence  analysis  for  smooth  solutions  and  in  steady  state  convergence. 
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4.  The  approximation  is  self  similar.  That  is,  when  fully  discrete  with  Runge-Kutta  methods  in  section 
2,  the  scheme  is  invariant  when  the  spatial  and  time  variables  are  scaled  by  the  same  factor. 

We  then  indicate  how  the  scheme  is  generalized  in  a  more  complex  situation: 

1.  For  scalar  equations  without  the  property  f'(u)  >  0,  one  could  use  a  flux  splitting 


f(u)=f+(u)+f  (u), 


df+(u ) 


>  0, 


(If  (u) 


<0, 


du  du 

and  apply  the  procedure  above  to  f+(u ),  and  a  mirror  image  (with  respect  to  j  +  1/2)  procedure 
to  /  (</).  The  only  requirement  for  the  splitting  is  that  /±(u)  should  be  as  smooth  functions  of 
u  as  /(«)  is  and  as  the  order  of  the  scheme  requires  (e.g.  if  the  scheme  is  fifth  order,  f(u)  and 
/±(u)  should  all  have  five  continuous  derivatives  with  respect  to  u).  In  most  applications  the  simple 
I  .ax  Friedrichs  flux  splitting 


1 


/±(m)  =  x(/(u)  ±  w),  a  =  maxu\f(u)\ 


where  the  maximum  is  taken  over  the  relevant  range  of  u,  is  a  good  choice. 

2.  For  systems  of  hyperbolic  conservation  laws,  the  nonlinear  part  of  the  WENO  procedure  (i.e.  the 

determination  of  the  smoothness  indicators  0k  and  hence  the  nonlinear  weights  w  -,)  should  be  carried 
out  in  local  characteristic  fields.  Thus  one  would  first  find  an  average  Uj+ 1/2  of  uj  and  uj+ 1  (e.g. 
the  Roe  average  [37]  which  exists  for  many  physical  systems),  and  compute  the  left  and  right 
eigenvectors  of  the  Jacobian  /'(uJ+1/2)  and  put  them  into  the  rows  of  RJ^^  and  the  columns 
of  Rj+ i/o,  respectively,  such  that  RJ^^  /'(wy+1/2)  =  A j+1/2  where  A j+1/2  is  a  diagonal 

matrix  containing  the  real  eigenvalues  of  f'(Uj+ i/o).  One  then  transforms  all  the  quantities  needed 
for  evaluating  the  numerical  flux  fj+1/2  to  the  local  characteristic  fields  by  left  multiplying  them 
with  R“)11^0,  and  then  computes  the  numerical  fluxes  by  the  scalar  procedure  in  each  characteristic 
field.  Finally,  the  flux  in  the  original  physical  space  is  obtained  by  left  multiplying  the  numerical 
flux  obtained  in  the  local  characteristic  fields  with  Rj+i  / 2- 

3.  If  one  has  a  non-uniform  but  smooth  mesh,  for  example  x  =  x(£)  where  is  uniform  and  x(£) 
has  at  least  five  continuous  derivatives  for  the  fifth  order  method,  then  one  could  use  the  chain 
rule  f(u)x  =  f(u)^/x'(^)  and  simply  use  the  procedure  above  for  uniform  meshes  to  approximate 
f(u ■)£.  The  metric  derivative  #'(£)  should  be  either  obtained  through  an  analytical  formula  (if  the 
transformation  x  =  x(£)  is  explicitly  given)  or  by  a  finite  difference  approximation  which  is  at  least 
fifth  order  accurate,  for  example  again  by  a  WENO  approximation.  Using  this,  one  could  use  finite 
difference  WENO  schemes  on  smooth  curvilinear  coordinates  in  any  space  dimension. 

4.  WENO  finite  difference  schemes  are  available  for  all  odd  orders,  see  [23]  and  [2]  for  the  formulas  of 
the  third  order  and  seventh  through  eleventh  order  WENO  schemes. 

We  present  two  numerical  examples  to  illustrate  the  capability  of  the  finite  difference  WENO  schemes. 
Both  are  obtained  with  the  fifth  order  WENO  schemes. 

The  first  example  is  the  double  Mach  reflection  problem,  originally  given  in  [45]  and  later  used  often  in 
the  literature  as  a  benchmark.  The  computational  domain  is  [0, 4]  x  [0, 1],  although  typically  only  the  results 
in  [0,3]  x  [0, 1]  are  shown  in  the  figures.  The  reflecting  wall  lies  at  the  bottom,  starting  from  x  =  Initially 
a  right-moving  Mach  10  shock  is  positioned  at  #  =  \,y  =  0  and  makes  a  60°  angle  with  the  x  axis.  For  the 
bottom  boundary,  the  exact  post-shock  condition  is  imposed  for  the  part  from  x  =  0  to  x  =  |  and  a  reflective 
boundary  condition  is  used  for  the  rest.  At  the  top  boundary,  the  flow  values  are  set  to  describe  the  exact 
motion  of  a  Mach  10  shock.  The  computation  is  carried  out  to  t  =  0.2.  At  a  very  refined  resolution,  the 
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Fig.  3.1.  Density  contours,  double  Mach  reflection,  fifth  order  finite  difference  WENO  scheme.  Left:  uniform  mesh  with 
Ax  =  Ay  =  Right:  non-uniform  moving  mesh  with  1/4  as  many  2D  mesh  points  (480  points  in  y). 


1921x481 
1=  0.2 


Fig.  3.2.  A  ‘‘zoomed  in"  version  of  the  density  contours,  double  Mach  reflection,  fifth  order  finite  difference  WENO 
scheme.  Left:  uniform  mesh  with  Ax  =  Ay  =  ;  Right:  non-uniform  moving  mesh  with  1/4  as  many  2D  mesh  points  (480 

points  in  y ) . 

slip  line  induced  instability  and  roll-up  can  be  observed,  see,  e.g.  the  adaptive  mesh  simulation  in  [3].  The 
capability  of  a  numerical  method  to  simulate  these  roll-ups  is  an  indication  of  its  small  numerical  viscosity 
and  high  resolution.  In  Fig.  3.1,  left,  we  give  the  density  contours  of  the  simulation  result  with  the  fifth 
order  WENO  scheme  on  a  fixed,  uniform  mesh  with  Ax  =  Ay  =  g^.  In  Fig.  3.1,  right,  we  give  the  density 
contours  of  the  simulation  result  with  the  fifth  order  WENO  scheme  on  a  non-uniform  and  moving  mesh, 
which  is  smooth  and  concentrates  its  points  near  the  shock  and  the  region  under  the  double  Mach  stem,  with 
only  half  the  number  of  points  in  each  direction  (480  points  in  y).  The  mesh  movements  were  determined 
by  a  given  smooth  transformation  which  follows  the  structure  of  the  solution.  Fig.  3.2  gives  a  “zoomed  in” 
picture.  We  can  clearly  see  that  the  resolutions  are  comparable  while  the  moving  non-uniform  mesh  version 
uses  only  1/4  as  many  2D  mesh  points  as  the  uniform  one,  hence  saving  a  lot  of  computational  effort. 

The  second  example  is  the  problem  of  a  supersonic  flow  past  a  cylinder  [23].  In  the  physical  space,  a 
cylinder  of  unit  radius  is  positioned  at  the  origin  on  the  x-y  plane.  The  computational  domain  is  chosen  to 
be  [0,1]  x  [0, 1]  on  the  £-/?  plane.  The  mapping  between  the  computational  domain  and  the  physical  domain 
is: 


*  =  (Rx  -  (Rx  -  m  C0S(£(2,?  -  1)),  y  =  (Ry  -  (Ry  -  1)£)  sin(0(2,?  -  1)), 

where  Rx  =  3 ,Ry  =  6  and  f)  =  jf.  Fifth  order  finite  difference  WENO  and  a  uniform  mesh  of  60  x  80 
points  in  the  computational  domain  are  used.  The  problem  is  initialized  by  a  Mach  3  shock  moving  towards 
the  cylinder  from  the  left.  Reflective  boundary  condition  is  imposed  at  the  surface  of  the  cylinder,  i.e.  at 
£  =  1,  inflow  boundary  condition  is  applied  at  £  =  0  and  outflow  boundary  condition  is  applied  at  =  0, 1. 
We  present  an  illustration  of  the  mesh  in  the  physical  space  (drawn  every  other  grid  line),  and  the  pressure 
contour,  in  Fig.  3.3.  We  can  clearly  see  that  the  finite  difference  WENO  scheme  can  handle  such  curvilinear 
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Fig.  3.3.  Flow  past  a  cylinder.  Left:  physical  grid  drawn  on  every  other  grid  line;  right:  pressure  contours  with  20  contour 
lines. 

meshes  very  well. 

In  summary,  we  can  say  the  following  about  finite  difference  WENO  schemes: 

1.  They  can  only  be  used  for  regular  geometry  that  can  be  covered  either  by  uniform  or  smooth 
curvilinear  meshes.  The  smoothness  of  the  mesh  must  be  comparable  with  the  order  of  accuracy  of 
the  scheme  in  order  to  obtain  a  truly  high  order  result. 

2.  If  the  computational  problem  allows  for  such  meshes,  then  the  finite  difference  WENO  schemes  are 
good  choices  as  they  are  easy  to  code  and  fast  to  compute,  especially  for  multi  dimensional  problems. 
Usually  the  fifth  order  WENO  scheme  is  the  best  choice,  unless  the  nature  of  the  problem  asks  for 
higher  orders  of  resolution. 
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3.  The  finite  difference  WENO  schemes  can  also  be  used  in  an  adaptive  mesh  environment,  provided 
that  a  smooth  (in  space  and  time)  mesh  can  be  generated.  To  generate  such  meshes  is  not  easy, 
especially  for  higher  order  schemes  where  the  requirement  for  the  smoothness  of  the  meshes  is 
stronger. 


4.  Finite  volume  WENO  schemes.  A  finite  volume  scheme  for  a  conservation  law  such  as  (1.1.1) 
approximates  an  integral  version  of  it.  Thus,  the  computational  domain  is  partitioned  into  a  collection  of 
cells  A i,  which  in  2D  could  be  rectangles,  triangles,  etc.,  and  the  cell  averages  of  the  solution  u 


Ui(t) 


1 

WJ 


u(x,  y,  t)dxdy 


are  the  numerically  approximated  quantities.  If  A j  is  a  control  volume,  the  semi-discrete  finite  volume 
scheme  of  equation  (1.1.1)  is: 


d  ,  ,  1 

5A(() + ra 


F  ■  n  ds  =  0 


sa4 


(4.4.1) 


where  F  =  ( f,g ),  and  n  is  the  outward  unit  normal  of  the  cell  boundary  dAj.  The  line  integral  in  (4.4.1)  is 
typically  discretized  by  a  Gaussian  quadrature  of  sufficiently  high  order  of  accuracy, 


F  ■  n  ds 


Q 

I5AJ'I  u kF(u(Gk,t ))  •  n 

k=l 


and  F(u(Gk,t))  ■  n  is  replaced  by  a  numerical  flux  (approximate  or  exact  Riemann  solvers).  For  example, 
one  could  use  the  simple  Lax-Friedrichs  flux,  which  is  given  by 

F(u(Gk,t ))  •  n  ^  [( F(u~(Gk,t ))  +  F(u+(Gk,t)))  ■ n-a  ( u+(Gk,t )  -  u~(Gk,t))] 

where  a  is  taken  as  an  upper  bound  for  the  eigenvalues  of  the  Jacobian  in  the  n  direction,  and  u~  and  u+ 
are  the  values  of  u  inside  the  cell  Aj  and  outside  the  cell  A  j  (inside  the  neighboring  cell)  at  the  Gaussian 
point  Gk. 

Clearly,  the  success  of  the  finite  volume  scheme  depends  crucially  on  a  good  “reconstruction”  procedure, 
which  is  the  procedure  to  obtain  high  order  and  non-oscillatory  approximations  to  the  solution  u  at  the 
Gaussian  points  along  the  cell  boundary,  u±{Gk,t),  from  the  neighboring  cell  averages.  Usually,  this  recon¬ 
struction  problem  is  handled  in  the  following  way:  given  a  stencil  of  R  =  cells,  find  a  polynomial  of 

degree  r,  whose  cell  average  in  each  cell  within  the  stencil  agrees  with  the  given  cell  average  of  u  in  that  cell. 
This  gives  a  linear  system  of  R  equations  and  R  unknowns  (the  coefficients  of  the  polynomial  when  expanded 
in  a  certain  basis),  and,  if  it  has  a  unique  solution,  the  polynomial  can  be  evaluated  at  the  Gaussian  point 
to  get  the  approximation  to  u±(Gk,t).  In  practice,  there  are  a  lot  of  complications  in  this  procedure,  as  not 
all  stencils  result  in  a  solvable  or  well  conditioned  linear  system.  One  would  often  resort  to  a  least  square 
procedure  with  more  than  the  necessary  number  of  cells  in  the  stencil  to  solve  this  problem,  see,  e.g.  [20].  If 
the  cells  are  rectangles  rather  than  triangles,  then  a  tensor  product  polynomial  and  a  tensor  product  stencil 
would  be  much  easier  to  work  with  [38]. 

A  typical  WENO  finite  volume  scheme  is  constructed  as  follows: 

1.  We  identify  several  stencils  Sj,  i  =  1 ..... r/.  such  that  the  control  volume  Aj  belongs  to  each  stencil. 

Q 

We  denote  by  T  =  (J  the  larger  stencil  which  contains  all  the  cells  from  the  q  stencils. 

2  —  1 
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2.  We  obtain  a  (relatively)  lower  order  reconstruction  polynomial,  denoted  by  pi(x),  associated  with 

each  of  the  stencils  Sj,  for  i  =  We  also  obtain  a  (relatively)  higher  order  reconstruction 

polynomial,  denoted  by  P{x ),  associated  with  the  larger  stencil  T- 

3.  We  find  the  combination  coefficients,  also  called  linear  weights,  denoted  by  71,  ...  ,  7g,  such  that  for 
a  Gaussian  point  Gk  on  the  cell  boundary, 

Q 

P(Gk )  =  ^7iPi(G*) 

i- 1 

for  all  possible  given  cell  averages  in  the  stencil.  These  linear  weights  depend  on  the  mesh  geometry, 
the  point  Gk,  and  the  specific  reconstruction  requirements,  but  not  on  the  given  cell  averages  in  the 
stencil. 

4.  We  compute  the  smoothness  indicator,  denoted  by  Si,  for  each  stencil  Si,  which  measures  how 
smooth  the  function  pi(x)  is  in  the  target  cell  A j.  The  smaller  this  smoothness  indicator  Si,  the 
smoother  the  function  pi(x)  is  in  the  target  cell.  These  smoothness  indicators  are  obtained  with  the 
same  integral  formulas  as  in  the  finite  difference  WENO  schemes.  The  details  can  be  found  in  [20] 
and  [38]. 

5.  We  compute  the  nonlinear  weights  based  on  the  smoothness  indicators: 


Wi  = 


Wj 

Efc  wk' 


(£  +  Sk)2 


where  7*  are  the  linear  weights  determined  in  step  3  above,  and  e  is  again  a  small  number  to  avoid 
the  denominator  to  become  0  and  is  usually  taken  as  e  =  10-6  in  the  computations.  The  final 
WENO  reconstruction  is  then  given  by 


Q 

u~(Gk)  =  ^2 wiPi(Gk )• 

i= 1 


We  summarize  the  properties  of  this  WENO  finite  volume  scheme.  For  more  details,  see  [20]  and  [38]. 

1.  For  2D  triangulation  with  arbitrary  triangles,  third  and  fourth  order  finite  volume  WENO  schemes 
are  available,  [20],  [38].  The  third  order  scheme  is  quite  robust.  The  fourth  order  scheme,  however, 
seems  to  have  more  restrictive  requirements  on  the  triangulation  for  stability  for  solving  systems  of 
conservation  laws. 

2.  For  2D  triangulation  with  tensor  product  rectangle  meshes,  which  could  be  non-uniform  and  non¬ 
smooth,  the  fifth  order  WENO  scheme  in  [38]  is  quite  robust  and  gives  very  good  numerical  results. 

We  will  again  use  the  double  Mach  reflection  problem  to  illustrate  the  behavior  of  the  finite  volume 
WENO  schemes.  To  save  space  we  will  show  only  the  results  obtained  with  the  fifth  order  finite  volume 
WENO  scheme  on  a  tensor  product  mesh,  with  a  uniform  mesh  of  A;r  =  Ay  =  ^  [38],  in  Fig.  4.1.  Results 
obtained  with  the  third  and  fourth  order  WENO  schemes  on  triangular  meshes  can  be  found  in  [20]. 

In  summary,  we  can  say  the  following  about  finite  volume  WENO  schemes: 

1.  They  can  be  used  for  arbitrary  triangulation.  However  they  are  much  more  complex  to  code  and 
much  more  expensive  in  CPU  cost  than  finite  difference  WENO  schemes  of  the  same  order  of 
accuracy.  This  is  because  they  have  to  rely  on  multidimensional  reconstructions  (polynomials  of  2 
or  3  variables  in  2D  or  3D) ,  and  the  flux  integrals  on  the  cell  boundaries  must  be  performed  by  multi 
point  Gaussian  quadratures.  As  a  rule  of  thumb,  a  finite  volume  WENO  scheme  is  at  least  4  times 
more  expensive  in  2D  and  9  times  more  expensive  in  3D,  compared  with  a  finite  difference  WENO 
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Fig.  4.1.  Double  Mach  reflection,  fifth  order  finite  volume  WENO  scheme,  uniform  mesh  with  Ax  =  Ay  =  Left: 

density  contours;  Right:  a  ‘‘zoomed  in”  version  of  the  density  contours. 

scheme  on  the  same  mesh  and  of  the  same  order  of  accuracy,  see,  e.g.  [5]  for  such  a  comparison  for 
ENO  schemes. 

2.  Finite  volume  WENO  schemes  on  a  tensor  product  mesh  are  more  robust  and  can  be  constructed 
for  higher  order  of  accuracy  than  finite  volume  WENO  schemes  on  arbitrary  triangulation. 

3.  Finite  volume  WENO  schemes  should  be  used  in  the  situation  when  it  is  impossible  to  apply  a 
smooth  curvilinear  mesh. 

5.  Discontinuous  Galerkin  methods.  Similar  to  a  finite  volume  scheme,  a  discontinuous  Galerkin 
(DG)  method  for  a  conservation  law  such  as  (1.1.1)  also  approximates  an  integral  version  of  it.  The  compu¬ 
tational  domain  is  again  partitioned  into  a  collection  of  cells  A,,  which  in  2D  could  be  rectangles,  triangles, 
etc.,  and  the  numerical  solution  is  a  polynomial  of  degree  r  in  each  cell  At.  The  degree  r  could  change 
with  the  cell,  and  there  is  no  continuity  requirement  of  the  two  polynomials  along  an  interface  of  two  cells. 
Thus,  instead  of  only  one  degree  of  freedom  per  cell  in  a  finite  volume  scheme,  namely  the  cell  average  of 
the  solution,  there  are  R  =  ('-+i)b-+2)  jegrees  0f  freedom  per  cell  for  a  DG  method  using  piecewise  r-th 
degree  polynomials  in  2D.  These  R  degrees  of  freedom  are  chosen  as  the  coefficients  of  the  polynomial  when 
expanded  in  a  local  basis.  One  could  use  a  locally  orthogonal  basis  to  simplify  the  computation,  but  this  is 
not  essential. 

The  DG  method  is  obtained  by  multiplying  (1.1.1)  by  a  test  function  v(x,y)  (which  is  also  a  polynomial 
of  degree  r  in  the  cell),  integrating  over  the  cell  A j,  and  integrating  by  parts: 

—  I  ti(r,  I/.  I)r(.r.  yjd.rdy  j  F(u)  ■  Vu  dxdy  +  f  F(u)-nvds  =  0 

dt  J Aj  J Aj  JdAj 

where  the  notation  and  the  treatment  of  the  line  integral  are  the  same  as  in  the  finite  volume  scheme  (4.4.1). 
The  extra  volume  integral  term  fA  F(u)  ■  Vvdxdy  can  be  computed  either  by  a  numerical  quadrature  or 
by  a  quadrature  free  implementation  [1]  for  special  systems  such  as  the  Euler  equations  (1.1.1).  Notice  that 
if  a  locally  orthogonal  basis  is  chosen,  the  time  derivative  term  JA  u(x,  y,t)v(x,  y)dxdy  would  be  explicit 
and  there  is  no  mass  matrix  to  invert.  However,  even  if  the  local  basis  is  not  orthogonal,  one  still  only  needs 
to  invert  a  small  R  x  R  local  mass  matrix  (by  hand)  and  there  is  never  a  global  mass  matrix  to  invert  as  in 
a  typical  finite  element  method. 

When  applied  to  problems  with  smooth  solutions,  the  DG  method,  as  briefly  described  above,  can 
already  be  used  as  is.  For  problems  containing  discontinuous  solutions,  however,  a  nonlinear  total  variation 
bounded  (TVB)  limiter  might  be  needed.  For  details,  see  [39,  6,  8,  9]. 

We  summarize  the  properties  of  the  DG  method  here.  For  more  details,  see  [11]. 
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Rectangles  P2,  Ax  =  Ay  =  1/480 


Fig.  5.1.  Double  Mach  reflection,  third  order  ( P 2  polynomials)  discontinuous  Galerkin  method,  rectangular  uniform  mesh 
with  Ax  =  Ay  —  .  Left:  density  contours;  Right:  a  ‘‘zoomed  in”  version  of  the  density  contours. 

1.  The  DG  method  has  the  best  provable  stability  property  among  all  three  methods  discussed  in  this 
paper.  One  can  prove  a  cell  entropy  inequality  for  the  square  entropy  [22],  which  implies  L2  stability 
for  the  full  nonlinear  case  with  possible  discontinuous  solutions,  and  any  converged  solution  is  an 
entropy  solution  for  a  convex  scalar  conservation  law.  This  cell  entropy  inequality  holds  for  all  scalar 
nonlinear  conservation  laws,  all  orders  of  accuracy  of  the  scheme,  all  space  dimensions,  arbitrary 
triangulation,  and  without  the  need  to  use  the  nonlinear  limiters. 

2.  The  DG  method  can  also  be  used  on  problems  with  second  derivatives  (diffusion  terms  such  as  those 
from  the  Navier-Stokes  equations),  [10],  [34].  It  can  even  be  used  on  problems  with  third  derivative 
terms  [46].  Theoretical  results  about  stability  and  rate  of  convergence  are  very  similar  to  those  for 
the  first  derivative  PDEs.  Unlike  the  traditional  mixed  method,  such  local  discontinuous  Galerkin 
methods  for  higher  derivatives  are  truly  local  (the  auxiliary  variables  introduced  for  the  derivatives 
can  be  eliminated  locally)  and  share  with  the  discontinuous  Galerkin  method  all  the  flexibility  and 
advantages  such  as  a  tolerance  of  arbitrary  triangulation  with  hanging  nodes,  parallel  efficiency, 
easiness  in  h-p  adaptivity,  etc. 

We  will  again  use  the  double  Mach  reflection  problem  to  illustrate  the  behavior  of  the  DG  methods.  We 
present  the  result  of  the  third  order  method  (piecewise  quadratic  polynomials)  on  a  rectangular  mesh  with 
Ax  =  Ay  =  Aq  [9],  in  Fig.  5.1. 

In  summary,  we  can  say  the  following  about  the  discontinuous  Galerkin  methods: 

1.  They  can  be  used  for  arbitrary  triangulation,  including  those  with  hanging  nodes.  Moreover,  the  de¬ 
gree  of  the  polynomial,  hence  the  order  of  accuracy,  in  each  cell  can  be  independently  decided.  Thus 
the  method  is  ideally  suited  for  h-p  (mesh  size  and  order  of  accuracy)  refinements  and  adaptivity. 

2.  The  methods  have  excellent  parallel  efficiency.  Even  with  space  time  adaptivity  and  load  balancing 
the  parallel  efficiency  can  still  be  over  80%. 

3.  They  should  be  the  methods  of  choice  if  geometry  is  complicated  or  if  adaptivity  is  important, 
especially  for  problems  with  smooth  solutions. 

4.  For  problems  containing  strong  shocks,  the  nonlinear  limiters  are  still  less  robust  than  the  advanced 
WENO  philosophy.  There  is  a  parameter  (the  TVB  constant)  for  the  user  to  tune  for  each  problem. 
For  rectangular  meshes  the  limiters  work  better  than  for  triangular  ones.  Other  limiters  are  still 
being  investigated  in  the  literature. 

6.  Concluding  remarks.  We  have  discussed  three  classes  of  typical  high  order  numerical  methods 
used  in  CFD,  especially  for  problems  containing  both  shocks  or  high  gradient  regions  and  complex  smooth 
region  structures.  These  are  finite  difference  WENO  schemes,  finite  volume  WENO  schemes  and  discontin- 
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uous  Galerkin  methods.  All  three  methods  use  the  same  nonlinearly  stable  high  order  Runge-Kutta  time 
discretizations  [43],  hence  their  difference  is  only  in  spatial  discretizations.  Finite  difference  WENO  schemes 
have  the  advantage  of  simplicity  and  lower  CPU  cost,  especially  for  multi  dimensional  problems,  but  they 
can  only  be  applied  on  smooth  structured  curvilinear  meshes.  If  a  computational  problem  allows  the  usage 
of  such  meshes,  finite  difference  WENO  schemes  are  good  choices.  In  this  class  the  one  used  most  often  is 
the  fifth  order  WENO  scheme  in  [23].  Finite  volume  WENO  schemes  are  more  expensive  than  their  finite 
difference  counter  parts.  However,  they  do  have  the  advantage  of  allowing  arbitrary  triangulation,  at  least 
in  principle.  For  two  dimensional  triangulation  with  arbitrary  triangles,  WENO  finite  volume  schemes  of 
third  and  fourth  order  accuracy  are  available  [20],  [38].  The  third  order  version  is  quite  robust,  however  the 
fourth  order  version  seems  to  have  more  restrictive  requirements  on  the  type  of  triangulation  for  stability. 
Higher  order  versions  and  three  dimensional  cases  are  still  under  development.  For  structured  meshes,  finite 
volume  WENO  schemes  of  fifth  order  accuracy  [38]  are  very  robust  and  allow  for  arbitrary,  non-smooth  mesh 
sizes,  hence  they  can  be  used  in  more  general  situations  than  the  finite  difference  WENO  schemes.  Finally, 
the  discontinuous  Galerkin  method  is  the  most  flexible  in  terms  of  arbitrary  triangulation  and  boundary 
conditions.  It  is  ideally  suited  for  problems  with  smooth  solutions.  For  problems  containing  shocks,  the 
total  variation  bounded  limiter  [39,  6,  8,  9]  works  quite  well  for  rectangular  meshes,  and  reasonably  well  for 
arbitrary  triangulation.  However  they  are  still  less  robust  than  WENO  schemes  as  they  contain  a  tuning 
parameter.  An  active  research  direction  now  is  the  search  for  a  more  robust  and  high  order  preserving  limiter 
for  general  triangulation. 
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